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Abstract. We study the finite-size effect in calculating the autocorrelation function of the 
global flux using periodic boundary conditions in this Letter. By studying the spatiotemporal 
correlation of the local flux, it is shown that the size effect is induced by the interaction of the 
sound modes. A local-flux fluctuation may excite sound modes spreading outwards with a speed 
v, and interaction may take place when they meet due to the finite size of simulation domain. As 
a consequence, we conclude that one could only correctly extract the correlation function using a 
periodic domain of L for a time period of t < L/(2v). This criterion implies that a lot of previous 
works related to the numerical simulation of the correlation function need to be rechecked. 
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1. Introduction 

The correlation functions of fluctuations are essential in exploring transport behaviors. Particularly, 
applying the auto-correlation functions of flux fluctuations in equilibrium systems, one can extract 
the corresponding transport coefficients using the Green-Kubo formula For a one-dimensional 

(ID) system, for example, one can apply 

1 f 

K = lim lim / Cjj(t)dt (1) 

T-»ooL->oo Kgi^L Jq 

to calculate the thermal conductivity k [Bl-[l4l|. where Cjj(t) = (J(O)J(i)) is the auto-correlation 
function of fluctuations of the global heat flux J(t). In the formula, r is the time of evolution; L is 
the system size; ks is the Boltzmann constant; T is the temperature of the system. Predicting the 
thermal conductivity using this formula becomes practically important in recent years because of the 
appearance of the nanoscale materials and devices while the experimental measurement of thermal 
conductivity for these low-dimensional systems is quite difficult. For this purpose, the asymptotic 
behavior of the correlation function is particular essential since it determines the convergence of 
the conductivity. When the correlation function decays in a power-law way, the index of power 



determines further the divergent law of the conductivity on the system size [lCj, lll| . 

In principle, the correlation function should be extracted with an infinite L. In practical 
simulations, the size of a system is finite. The conventional strategy for this situation demands that 
the size should be large enough. However, there is still no clarification how large is enough to avoid 
the finite-size effect. 

In applying the fixed and open boundary conditions, the correlation function is found oscillating 
as a function of time [l(| . The phenomenon is considered to be finite-size effect due to the reflection 
of the boundaries. These kinds of boundaries are thus unsuitable for extracting the correlation 
function. In applying the periodic boundary condition, the system is equivalent to a circle with 
length L. The flux flows along the circle and the boundary reflection is avoided. Periodic boundaries 
thus effectively overcome the limitation of a small simulation domain and usually adopted in 
calculating the correlation function of the flux Q. However, oscillation phenomenon with smaller 
amplitude is also observed in gas models 12- T^. In lattice models, though oscillation phenomenon 



has not been observed, it is found that the accuracy of thermal conductivity is sensitive to L 0-0| • 
Meanwhile, direct simulations indicate that size effect do appear for either the gas or lattice models 
since the correlation function obtained with smaller periodic domain drops from that obtained with 
large-size domain after a certain period of time. 

The open problems related to exploring the correlation function are as following. Why size 
effect can not be avoided by using periodic boundary condition? What is the underlying microscopic 
mechanism behind the oscillation phenomenon? Could the oscillation phenomenon induce the 
dissipation of the flux? In addition, whether one can explore the asymptotic behavior of the 
correlation function by studying the long-time behavior of a finite system as been done in previous 
references has not been checked. And a criterion for judging the degree of approximation 



of the correlation function obtained with a finite L to that with infinite size is still absent. 

In this paper we study two ID models, one is a gas model and the another is a lattice model, 
to clarify the above questions. We demonstrate firstly that the periodic oscillation phenomenon of 
the autocorrelation function of the global flux existing not only in the gas model but also in the 
lattice model with periodic boundary conditions. We then argue that the autocorrelation function 
of the global flux can be represented by the sum of the spatiotemporal correlation functions of the 
local flux. By studying the local flux, the mechanism of the periodic phenomenon is understood 
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Figure 1. Energy current auto-correlation function Cjj(t)/Cjj(0) versus t at periodic boundary 
condition, (a) corresponding to that of the gas model. The inset in (a) is a zoom of the same 
plot, (b) corresponding to that of the FPU model. The inset in (b) is the plot of the same data 
multiplied a rescaling factor t 0,67 . 



and a time criterion for approaching the autocorrelation function of the global flux using periodic 
boundary conditions is presented. 

2. Models 

The gas model (l2l . [20l | is consists of N hard-core point particles arranged in order in a ID box 
of length L = N with alternative mass mi for odd-numbered particles and ttt-2 for even-numbered 
particles. The particles travel freely except elastic collisions with their neighbors. Without loss of 
generality, we fix mi — 1, TO2 = 3 and set the energy density to be unity, as adopted by Cipriani et 
al [2(|. The lattice model is the so-called Fermi-Pasta- Ulam (FPU) model [2l[ with the Hamiltonian 

H = Y,(pW + V(x 1 -x^ 1 ) (2) 

i 

with V(x) — x 2 /2 + x 4 /A, where Xi is the displacement of the ith particle from its equilibrium 
position and pi is its momentum. In the simulations particles are assumed to have unit mass and 
the energy density is set to be unity, too. 

The global heat flux is identical with the global energy flux, since the total momentum is set 
to be zero and there is no residual global velocity [9j. Thus, the global heat flux in the gas model 
is calculated by J = J2 m i v i/^ For lattice models, we apply the conventional definition [22j 
J = Xi§j- to calculate the total heat flux, where the sum is over all the particles. 
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3. Results 

Figure 1 shows the normalized correlation function of the global-flux fluctuations for the two models 
with periodic boundary conditions. The oscillation is directly visible for the gas model, while 
periodic feature is invisible for the lattice model in the original scale. However, after a properly 
rescaling, as shown the inset in Figure 1(b), the periodic pulses appear, each pulse having an 
upwards peak followed by a downwards one. 

It has been shown that the spatiotemporal correlation function of the local-flux fluctuations 
can reveal important information of the correlation function of the global flux (3, El, 17-lj|. Here 



we demonstrate that it provides essential clues for understanding the periodic oscillation behavior 
on Cjj(t). The chain is divided into iVj, bins. Then the energy flux J can be expressed as the sum of 
the local flux, i.e., J = £jfc, where jk represents the local flux of the fcth bin. The spatiotemporal 

k 

correlation function, 

C(x,t) = {jk(0)ji(t)), (3) 

defines the correlation of the fluctuation jk(0) in the kth bin to the fluctuation of ji(t) in the Zth 
bin at time t. where x = (I — k)h , and h is the width of a bin. The time evolution of this function 
reveals the transport process of fluctuations of the local flux. 

Figure 2 shows the evolution of C(x,t) for the two models. In either cases, C(x,t) has the 
two-peak structure, indicating that the perturbations induced by the local-flux fluctuation diffuse 
outward from the initial bin. The induced fluxes are positively correlated to the initial flux since the 
peaks on C(x,t) keeps upwards. In another words, they have the same direction vector with jfc(O), 
though spread outward oppositely. The moving speed of peaks depends on the temperature (or 
energy density) of the system. In our simulations the energy density is unity for both models. We 
measured that v rs 1.75 for the gas model and v « 1.50 for the lattice model. They are independent 
on the amplitudes of the peaks. 

The two peaks should be sound modes predicted by the hydrodynamic theory 0,13]. The sound 
modes appear as two area-conserved peaks on the spatiotemporal correlation function of the local 
momentum fluctuations and energy fluctuations fl9l . ]23| , while as area-decaying peaks here since 
the momentum and energy are conserved quantities but the flux is not. The two-peak structure on 
the correlation function of local flux has been revealed in the gas model [13] . 

Let's take the gas model as an example to confirm that the peaks are sound modes. The sound 
velocity of an adiabatic gas is known to be v soun d = \J cpkpT / cym, where cp (cy) is the specific 
heat capacity at a constant pressure (volume) respectively and m is the molecular mass. In our gas 
model cp/cy — 3. Applying the average mass fh = (mi + ra^jl to be the molecular mass m and 
the Boltzmann constant ks = 1, we obtain v soun d ~ 1.73 in the case of unity energy density, which 
is consist to the speed measured on C(x, t) within the simulation error. 

Since (J(O)J(t)) = (£ Jfc(0) £ jl(t)) and the term on the right-hand side can be written as 

(Ejfc(O)EiiW) = E(ifc(0)£i*(*)>. we obtain 
fc 

{J(Q)J(t)) = ^Ck{x,t)das t (4) 
fc 

where the sum is over the bins. Each bin is statistically identical in the case of periodic boundary 
conditions. Therefore, the time-dependent behavior of the area of C(x, t) represents that of the 
correlation function Cjj(t) of the global flux. 
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Figure 2. The evolution of C(x,t), (a) for the gas model, (b) for the FPU model. 

The two peaks move freely before they meet at the opposite position for t < L/2v. During this 
period, Cjj(t) should behave as in the case of infinite system. At t — nL/2v with n = 1,2, 3. ..the 
two peaks meet again and again. When n is odd they meet at the opposite position to where they 
are excited, and when k is even they meet at the original position. If there is no interaction when 
they meet, they will pass across and result no effect on Cjj(t), and Cjj(t) should be the same as 
in the infinite system. 

To check whether they interact when they meet, Figure 3 shows the area of C(x,t) as a 
function of time. If there is no interaction, the area should decay monotonously. It can be seen 
that, nevertheless, in the gas model, the area drops around t = (2n — l)L/2i> while arises around 
t = 2nL/2v with n — 1,2,3,.... This fact indicates that interference takes place when two peaks 
cross through each other. As sound modes excited at the same source in an nonlinear system, it 
is not surprise that coherence may occur when they meet. However, we can not explain why they 
appear as constructive interference around the position that they are excited and as destructive 
interference around the opposite position. 

For the FPU model, the area of C(x,t) change almost monotonously with time. However, by 
proper;y rescaling, pulses with period L/2v appears. This fact implies that the mechanism of the 
interaction is different from that in the gas model. There is no global constructive and destructive 
interference in this model. 

The periodic behavior of C(x,t) matches well to that on Cjj(t). The periods are 294.3 on 
Cjj(t) and 295.0 on C(x,t) for the gas model, and are 44.0 on Cjj(t) and 44.0 on C(x,t) for the 
lattice model, respectively. Notice that the system sizes applied in the figures are L — 512 and 
L = 128 for these two models respectively, the periods can be well described by L/v = 292.6 and 
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Figure 3. (a) The time dependent behavior of the area of C(x, t) for the gas model under periodic 
boundary conduction; (b) the same but for the FPU lattice model, the insert is the same data 
multiplied a rescaling factor t o s exp(0.006t). 



L/2v = 42.7 correspondingly. Therefore, the periodic behavior on Cjj(t) should induced by the 
interaction of the sound modes. 

The essential question is then whether the repeating interaction of the sound modes would 
change the time-dependent behavior of Cjj(t). If the interaction induces no extra dissipation of 
the flux, the periodic oscillation or pules may exist but the time-dependent tendency of Cjj{t) 
should keep unchanged and one could still extract the decay law by smoothing the oscillations. 
Otherwise, Cjj(t) is reliable only before the collision of the peaks. To find the answer, we preform 
direct simulations. Figure 4 shows the dependence of the time-dependent behavior of Cj j(t) on 
the system sizes. This is a common way to study the size effect of the correlation f unct ions [ 1 (J-l 1 3| . 
Why and when a finite-size system fails to give the correct result has not been discussed. From the 
Figure 4(a), it can be seen that Cjjit) drops when the oscillation takes place, indicating that the 
interaction of the peaks in the gas model do arise the dissipation of the flux. In the case of the lattice 
model, Cjj(t) with smaller system size also drops from that with large sizes, seen Figure 4(b). For 
both models, the smaller the system size, the faster the decrease of Cjj(t). The phenomenon can 
be interpreted as the result of that the two peaks interact more frequently in smaller systems. As 
a consequence, Cjj(t) obtained by simulation is reliable only before the two peaks first meet at 
t = L/2v. 

The mechanism of the periodic behavior of the correlation function at the dynamics level should 
be topics in the future studies. Here we just briefly mention that if we accept the option that the 
sound modes in the FPU model is dominated by the solitary waves (l9l . |26| , the pulses on figure 
2 and figure 3 can be explained intuitively. Solitary waves can be excited in the FPU model is 
well-known (i3-[27|. Fi gure 5 shows such an example. A momentum kick with pi = 5 on a particle 



excite two solitary waves moving in opposite directions. When two solitary waves collides, besides 
the scattered solitary waves, a wave pocket is excited. The wave pocket is an effect of the collision 
of the solitary waves, indicating the extra dissipation to the energy and momentum as will as the 
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Figure 4. Heat current auto-correlation function Cjj(t)/N with various system sizes N for (a) 
the gas model and (b) the FPU lattice model. 



flux may arise in a finite system. Furthermore, it is known that when two solitary waves approaches 
close enough, they will accelerate while when they departs deceleration occurs [27J. This behavior 
may explain the pulse in the correlation of flux. 

4. Conclusions 

A fluctuation in an equilibrium can induce sound modes. The sound modes travel outwards in 
ID systems, appearing as two peaks on the spatiotemporal correlation functions of fluctuations 
of physical quantities such as momentum, energy and local flux. In an infinite system, they 
are traveling with constant speed v and never meet each other. In practical simulations using 
periodic boundary conditions, however, the simulation domain is finite and thus they will meet 
again and again. We show that when such two peaks meet each other, interaction would take place. 
The interaction is not only responsible for the periodic oscillation behavior of the autocorrelation 
function of the global flux but also for dissipation of the flux. As a result, the time-dependent 
behavior of the autocorrelation function of finite-size system will diverge from the function of an 
infinite-size system after the collision of the two peaks. Therefore, one can only extract reliable 
correlation function of flux by simulation using periodic boundary conditions in the time period of 
t < L/2v. 

As a consequence, the asymptotic behavior of the correlation function can not be approached 
with a finite-size system by extending the evolution time. This reminder is not be dispensable since 
one can find frequently in references that some researchers try to extract the asymptotic behavior 
in this way |28j|. In these studies, one may employ, say a finite system with L = 20000 to extract the 
time-dependent behavior of Cjj(t) up to t = 10 6 . At this time scale, repeatedly interaction of the 
sound modes has taken place and the decaying behavior of finite-size system should divert already 
from that of infinite-size system. The asymptotic behavior is essential to check the hydrodynamical 
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Figure 5. The collision process of two solitary waves of FPU model. Solitary waves excited by a 
kink with pi = 5 on the center particle of the FPU model. The other particles keep stationary at 
their equilibrium positions in the beginning, (a) a snapshot at t=400; (b) a snapshot at t=800. 



or kinetic predictions of the size-dependent law of the heat conductivity [id . and thus should 
be as precise as possible. 

We would like to point out that in calculating the heat conductivity using the Green-kubo 
formula researchers also apply a criterion of r = L/v to evolute the integral (1). This is a technique 
treatment in order to obtain the size-dependent behavior of the heat conductivity on the system 
size. Here we discuss the problem how to obtain the asymptotic behavior of the correlation function 
using the finite-size simulation. We show that only the function obtained in a time-period with 
t < L/2v is consistent with that of an infinite system. To approach the asymptotic behavior, one 
should study systems with sufficient size and investigate the behavior within this period. Indeed, 
the cutoff at r = L/v is not correct according our criterion, since the peaks have already experienced 
two times of collision during this period. 

We have not discuss the exact reliable time period since it is system-dependent. It can be 
defined by t = (L — l)/2v, where I is the width of the peaks on C(x,t) at t — L/2v. Before 
t = (L — l)/2v the two peaks have not meet at all and thus interaction has not involved. The width 
of the peaks and their time-dependent behavior depends on specific models. The mechanism of the 
dispersion behavior of the peaks should be understood in the future at the microscopic dynamics 
level. 
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